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We model the unsteady evolution of turbulent buoyant plumes following temporal changes 
to the source conditions. The integral model is derived from radial integration of the gov¬ 
erning equations expressing the conservation of mass, axial momentum and buoyancy in 
the plume. The non-uniform radial profiles of the axial velocity and density deficit in the 
plume are explicitly described by shape factors in the integral equations; the commonly- 
assumed top-hat profiles lead to shape factors equal to unity. The resultant model for 
unsteady plumes is hyperbolic when the momentum shape factor, determined from the 
radial profile of the mean axial velocity in the plume, differs from unity. The solutions of 
the model when source conditions are maintained at constant values are shown to retain 
the form of the well-established steady plume solutions. We demonstrate through a linear 
stability analysis of these steady solutions that the inclusion of a momentum shape factor 
in the governing equations that differs from unity leads to a well-posed integral model. 
Therefore, our model does not exhibit the mathematical pathologies that appear in pre¬ 
viously proposed unsteady integral models of turbulent plumes. A stability threshold for 
the value of the shape factor is also identified, resulting in a range of its values where the 
amplitude of small perturbations to the steady solutions decay with distance from the 
source. The hyperbolic character of the system of equations allows the formation of dis¬ 
continuities in the fields describing the plume properties during the unsteady evolution, 
and we compute numerical solutions to illustrate the transient development of a plume 
following an abrupt change in the source conditions. The adjustment of the plume to 
the new source conditions occurs through the propagation of a pulse of fluid through the 
plume. The dynamics of this pulse are described by a similarity solution and, through 
the construction of this new similarity solution, we identify three regimes in which the 
evolution of the transient pulse following adjustment of the source qualitatively differ. 


1. Introduction 

Turbulent buoyant plumes occur in numerous industrial and envir onmental settings 
llWoods 2010). Industrial examples includ e ventilation and heating ( e.g. Baines &: Turner 
1969 : Linden. Lane-Serff fc Sme^ll990), industrial chimn eys (e.g. ISlawson fc Csanadv 
1967 1 and waste-water disposal (e.g. iKoh fc Broo^ 19751. In the natural envi ronment, 
turbulent plu mes are found in meteorological (e.g. lE manuel| lR94ISteven ^ l2nn,5h oceano- 
graphical (e .g. Sneer fc Rona Il989l: IStraneo fc Cenedesel 2015ll . and volcanological (e.g. 
Woodsl[r988t Sparks et aLlll997 : Woodhouse et al. 2013[l settings. Model s of steady plume s. 


based on the integral modelling a p proach pioneered bvIZeldovichl (1193711. iSchmidtl (1194111. 
Rouse. Yih fc Humphrevsl ( 1952ll . IPriestlev fc Balll ( 1955 1 and Morton. Tavlor fc Turner 
( 1956[ l. have been applied extensively to understand plume dynamics. In this approach, 
the turbulent flow in the plume is described on a time-scale that is longer than the eddy 
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turnover time (the time scale that characterises the turbulent motion), and therefore 
the complicated turbulent motions are not explicitly modelled and only the evolution of 
mean flow quantities with distance from the source are described. A further simplihcation 
is obtained by integrating the mean flow quantities over the cross section of the plume, 
resulting in a system of nonlinear ordinary differential equations that describe the spatial 
development of the steady integral flow quantities such as the fluxes of mass, momentum 
and buoyancy. 

To obtain the integral fluxes, it is necessary to specify the radial and azimuthal de¬ 
pendence of the three-dimensional mean fields. When the flow domain and ambient con¬ 
ditions do not introduce an asymmetry, the time-averaged flow quantities are swirl free 
(i.e. there is no azimuthal dependence of the mean flow quantities). Experiments suggest 
that, at sufficiently high Reynolds number and for an unstratified ambient fluid, buoy¬ 
ant plumes attain self-similar radial profiles for both the mean axial velocity and the 
concentration of the species that gives rise to the density di fference between the plume 
and the ambient fluid sufficiently far from t he fluid source ( Pa j3amcqlaou_fc List! 1988 ; 
Shabbir fc Georg j 1994 : Wang fc Law 2002 : Ezzamel. Salizzoni fc Hundl2015ll . Further¬ 
more, higher-order turbulent quantities such as the turbulent intensity (the root mean 
square of velocity component fluctuations and concentration fluctua tions from the mean) 
and turbulent stresses al s o evolve with self-sim ilar radial profiles (iPapanicolaou fc List 
1988t I Wang fc Lawll200^ lEzzamel et al. 2015[l . The cross-sectional integration then re¬ 
sults in a system of ordinary differential equations that describes the evolution of the 
mass flux, momentum flux and buoyancy flux with distance from the source. For stratified 
ambient environments, the detailed structure and evolution of the mean and turbulent 
quantities have been scrutinised less thoroughly. However, the applicability of integral 
models assuming similarity of the mean flow profiles has been de monstrated through com¬ 
paris on of model predictions to laboratory ex periments (see e.g. Morton_eta/. Ill95fil: Ibisd 
1982) and in fiel d-scale applications (see e.g. Turned 1986; WoodsI 19881: Speer fc Rona 
1989[lKavel[2nnil . 

While averaging over the turbulent timescale and integrating over the plume cross- 
section significantly simplifies the mathematical description of the plume dynamics, de- 
tailed information a bout the turbulent flow is not fully captured. The integral model of 
Morton et all ( 19561) incorporates the turbulent nature of the flow through a parameter¬ 
ization of turbulent mixing as an inflow of fluid from the ambient to the plume, referred 
to as entrainment. The velocity of the entraining flow is linearly related to the mean axial 
velocity scale of the plume (as dimensional analysis demands, since the only velocity scale 
that remains following time-averaging and cross-sectional integration is the mean axial 
velocity of the plume). [Morton et al. ( 1956h take a constant entrainment coefficient, and 
the application of this model has been successful in des cribing steady buoyant plumes 
over a wide range of scales, from the labora tory (see e.g. iMorton et i7i.lll956 : List 1982) 
to plumes from large volcanic eruptions (e.g. IWoodall988l) . illustrating the success of the 
integral modelling approach with a constant entrainment coefficient. 

Detailed examination of laboratory experiments of turbulent plumes suggests that the 
entrainment coefficient does not have a constant universal valu^ but rather evolves as 


the f low develops ( Wang fc Law 20021: Kaminski. Tait fc Carazza 20051: Ezzamel et al 


2015[) . In particular, for plumes that are strongly forced with a flux of momentum at 
the source (referred to as buoyant jets), the entrainment coefficient transitions from 
a value appropriate for jets to the value for plumes as the flow becomes increasingly 
driven by buoyancy. Integr al models that include an evolving entrainment coefficient 
have been proposed (see e.g. 


nave been proposed (see e.g. inoxlll97C ; [Kaminski et all 2005): ICarazzo. Kaminski fc Taf3 
20061: ICraske fc van Reeuwiik 2015oJ l6). building on the integral modelling approach of 









































































































Unsteady turbulent buoyant plumes 3 

Priestley fc Balll ( 1955h whereby an integral expression for the conservation of axial ki¬ 


netic energy is used with conservation of momentum to derive an integral expression 
for conservation of mass. (This modelling approach is discussed further in appendix 1X1 1 
However, the assumption of a constant entrainment coefficient capt ures the lead ing order 
behaviour of turbulent buoyant plumes over a wide range of scales ( TurneI^ll986^ ■ and for 
fully-developed turbulent plumes sufficiently far from the source, a constant entrainment 
coefficient is appropriate and represents the similarity of the flow profiles and the turbu¬ 
lent e ntrainment processes at different heights in the plumes ( Turneilll986 : Ezzamel et al 

20ii). 


The applicability of steady models to describe the inherently unsteady, turbulent mo¬ 
tion relies on a separation of time scales. If the conditions at the plume source and in the 
ambient are held steady, the steady integral models well describe the plume be haviour on 
time scales that are long compared to the eddy turn-over time ( Woodd 20ldl . However, 
if the source conditions change on a time scale that is longer than the turbulent fluc¬ 
tuations, then a signature of the source v ariation may be seen in time-averaged plum e 
dynamics downstream of the source (e.g. Scase. Caulfield fc Dalziell 20081: IScas 3 1200911 . 
Unsteady sources occur frequently in natural settings, for example as the strength of a 
volcanic source changes in magnitude during an eruption. In addition, temporal changes 
in ambient conditions on a time scale similar to the ascent time of a fluid parcel are likely 
to result in a transient response of the plume. 


To model unsteady plumes, integrals models have been proposed (IDelichatsiosI 1979: 

Yu 1990: Vul’fson & BorodinI 200ll: Scase et al. 20066: Scase. Caulfield & Dalziel 2006a: 

Scase et a/.ll200Sl: Scase. As 

pden & Caulfieldl 20091: 

Scasell2009ri that extend the modellins 

approach of iMorton et al\ ( 

1956) while retaining some of the underlving assumptions of 


the steady model. In particular, the unsteady models capture the variations of flow quan¬ 
tities on a time-scale that is longer than the eddy turn-over time, and assume that the 
radial profiles of the mean axial velocity and the concentration of the buoyancy gener¬ 
ating species remain in a self-similar form throughout the evolution. However, given the 
difficulty in obtaining robust experimental results for plumes with time-varying source 
conditions, the underlying assumptions have yet to be scrutinised in detail. Numerical 
simulation of the governing equations can be used as a surrogate for laboratory ex¬ 
periment and allow de tailed investigation s of the turbulent flow properties throughout 
the modelle d domain ( Jiang fc Luo I 2 OOOI : IPlourde et al. 2008t Craske fc van ReeuwiikI 
2013t 201Aa ). The phy sical basis of the unsteady models have been further questioned by 


Scase fc Hewittl (2012 ) i n an a nalys is the stability of the steady solutions of the models of 
DelichatsiosI ( 19791 .IYuI (1990 1. and IScase et al\ (2006^ to small harmonic perturbations 


at the source. Scase fc Hewitt ( 2012ll find that the perturbations grow as they propagate 
through the plume, and that the growth rate increases without bound as the frequency 
of the source oscillations is increased; ther efore the models are ill-posed, suggesting a 
critical physical process has been neglected ( Joseph fc Sant 199n[l . 

In an attempt to ‘regularize’ the unsteady plume models, IScase &: Hewitd (120121) in¬ 
troduce a phenomenological diffusive term into the equation for conservation of axial 
momentum to represent turbulent mixing processes; it curtails the unbounded growth 
of short wavelength perturbations thus leading to a well-posed model. The form of the 
diffusive term has been investigated bv ICraske fc van Reeuwiik ( 2015a l for momentum- 
driven (non-buoyant) jets using direct nu merical simulations; the numerical simulations 
suggest that the diffusive term proposed bv IScase fc Hewittl (12012^ does not describe well 
the transien t evolution of momentu m-driven jets. Here we find that the diffusive term in¬ 
troduced bv IScase fc Hewitt ( 2012 ) leads to a new pathology in the system of equations; 
























































































































4 M.J. Woodhouse, J.C. Phillips & A. J. Hogg 

the steady states of the ‘regularized’ model are spatially unstable and therefore cannot 
be realised. 

We show here that the ill-posedness in the unsteady models analyzed bv IScase fc Hewitt 
1 2 OI 2 II is due to the assumption of a top-hat profile for the mean axial velocity (i.e. the 
axial velocity at any height is assumed to be radially invariant within the plume, and 
zero outside of the plume) and therefore a failure to account for the non-uniform ra¬ 
dial profile of the mean axial velocity of the plume. Non-uniform radi al profiles for the 
axial veloci ty were examined in the early studies of steady plumes by iPriestlev &: Ball 
1 1955h and Morton et al. 1 1956ll . However, solutions of steady plume models have the 
same form whether top-hat or non-uniform radial profiles for the mean axial velocity 
are adopt ed, as dimensional analysis demands, with on ly changes to coefficients in the 
solutions ( Morton et al. 19561 : LindenI 20001: Kavd 2008l l and many subsequent analyses 
of steady plumes have adopted the top-hat formulation. 

When a top-hat velocity profile is adopted, the cross-sectionally averaged mass and 
momentum of the plume are transported with the same rate, but the transport rates 
differ when the radial profile of the mean axial velocity is non-uniform. We therefore 
propose an integral model of unsteady plumes that explicitly accounts for the differ¬ 
ent transport rates of the cross-sectionally averaged mass and axial momentum of the 
plume, by introducing a ‘shape factor’ in the equation for the conservation of momentum 
that differs from unity when non-uniform radial profiles of the mean axial velocity are 
assumed. A shape factor that differs from unity in shallow-water hydraulic models has 
been shown to fundamentally alter solutions of the system of equations due to a change 
in the characteristics of the hyperbolic system ( Hogg fc Pritchai^l2004l l. Here we show 
that including a shape factor changes the character of the system of equations describ¬ 
ing unsteady plumes and leads to a well-posed system of equations without the need to 
include diffusive terms. 

A similar approach to regularisin g an unsteady integral model o f momentum driven 
jets has been proposed recently bv ICraske fc van Beeuwiik ( 2015 qI b ). Th rough analy¬ 
sis of their direct numerical simulations. ICraske fc van Beeuwiik (201530) develop a 
well-posed integral model of unsteady jets that includes a description of dispersion in 
the plume, resulting from the non-uniform radial profile of the mean axial velocity that 
leads to different transport rates for mass, axial momentum a nd kinetic energy (re¬ 
ferred to as type I dispersion by ICraske fc van ReeuwiikI 2O15oJ 0H. and the deviation 
of the radial profile of mean axial velocity from a self-similar form (referred to as type 
H dispersion). Type H dispersion is important in jets, where the flow structure evolves 
significantl y in the neighbourhood of the source even for temporally in variant source 
conditions! Wang fc Law! 2002t Kaminski et al\ 2005t Ezzamel et al. 2015 1. but is likely 
less important for fu lly developed turbu lent plumes for which the radial p r ofiles a re 
in self-similar forms ( Ezzamel et al. 2015 1. However, Craske fc van Beeuwiik ( 2O15 qI 01 
show that explicitly accounting for the non-uniform radial profile of mean axial veloc¬ 
ity (type I dispersion) in the integral equations is critical to the well-posedness of the 
unsteady integral model for jets. 

In this contribution we demonstrate, through an analysis of the temporal evolution 
of small perturbations to steady solutions, that an integral model of unsteady plumes is 
well-posed when the momentum shape factor differs from unity. Furthermore, we identify 
a stability threshold in the value of the shape factor above which the amplitude of small 
perturbations decay. The system of equations we propose is hyperbolic, with a character¬ 
istic structure that, in certain situations, allows for the formation of ‘shocks’ during the 
transient evolution. Through the construction of similarity solutions, we identify scal¬ 
ing relationships that describe the propagation and growth of a transient pulse that is 
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advected through the plume following an abrupt change in the source conditions, and 
determine the regimes in which shocks are formed. 

This paper is organized as follows. In ^we provide a derivation of an integral model 
of unsteady buoyant plumes, beginning from the Reynolds-averaged Navier-Stokes equa¬ 
tions for a high Reynolds number flow together with an advection-diffusion equation for 
the transport of a scalar species the concentration of which determines the local fluid 
density. We demonstrate that the use of integral flow quantities, obtained through inte¬ 
gration over a plume cross-section, requires the inclusion of shape factors for the transport 
rates of the (cross-sectionally averaged) axial momentum and buoyancy of the plume. We 
show th at a momentum shap e factor modifies only slightly the classical power-law solu¬ 
tions of Morton et al\ (1 195611. In l}3l we re-examine the phenomenological diffusive term 
introduced bv IScase fc Hewitti (1^121 1 to ‘regularize’ the ill-posed unsteady plume model, 
and show this model introduces a new pathology into the system of equations whereby 
steady solutions are spatially unstable. We therefore examine the mathematical struc¬ 
ture and well-posedness of an unsteady integral model for plumes that accounts for the 
non-uniform radial profile of the mean axial velocity through a momentum shape factor, 
and we demonstrate in 21 that this leads to a well-posed system of equations. Numer¬ 
ical solutions of our unsteady model are presented in 21 to support the mathematical 
analysis. In 21 consider the evolution of a plume following an abrupt change in the 
source buoyancy flux and show that the adjustment of the plume occurs through the 
propagation of a pulse whose dynamics is described by a similarity solution. Finally, in 
21 we discuss the implications of our mathematical model and draw our conclusions. 


2. An integral model for nnsteady tnrbulent bnoyant plumes 

We model an unsteady turbulent buoyant plume formed due to the release of a fluid 
from a point source into an otherwise quiescent ambient fluid of a different density. A 
cylindrical coordinate system is adopted, with r and z denoting unit vectors in the 
radial and vertical directions, respectively. The plume and ambient are composed of 
incompressible fluids and the plume is assumed to have a circular (time-averaged) cross- 
section. The velocity field, m, is assumed to be axisymmetric, with u = ur wz. We 
further assume the plume is slender such that R/H 1 where R and H denote typical 
length scales in the radial and vertical direction respectively, that the Reynolds number 
of the emitted fluid at the source is sufficiently high such that the turbulent motion is 
fully developed throughout the flow domain, and that the plume fluid transports a scalar 
species (such as heat or salt) with the plume density linearly related to the concentration 
of the species. 

Turbulence in the plume is responsible for the entrainment of ambient fluid and mixing 
of the plume and ambient fluids, and its role is central to the ensuing dynamics. We 
therefore adopt the Reynolds-averaged Navier-Stokes equations to describe the fluid 
motion, with a Reynolds-averaged advection-diffusion equation to describe the transport 
of the scalar species that results in the density difference between the plume and the 
ambient fluids, taking u = u + u',w = w + w', and pr = ^ + where u and u' denote 
the ensemble average and fluctuation about the average of u, respectively, and similarly 
for w and pr, and where Pr = p (pa — p) /Po is the reduced gravity, with p denoting the 
gravitational acceleration, p and pa denoting the density of the plume and ambient fluids, 
respectively, and po is a characteristic density scale. The equations for the conservation 









6 M.J. Woodhouse, J.C. Phillips & A. J. Hogg 

of (the ensemble averaged) mass, axial momentum, and reduced gravity are then 


1 0 
r dr 


, X 9w 
(™) + -^ = 


dw 

dt 
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- — [ruw) + -^{w ) = gr 
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— — (ru'w') 
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dgr 

dt 


10, _ X 

+ - ^ [rgru) 
r dr 


dz 


(grw) 


Po 



dz 


{wpa) 


1 0 
r dr 


{rgr'u’) 


(2.1a) 


(2.16) 



(2.1c) 


Note, in (EH) we have made use of the plume slenderness and have invoked the Boussi- 
nesq approximation whereby differences in density are neglected except where they are 
multiplied by the gravitational acceleration. Furthermore, we have neglected molecular 
diffusion terms under the assumption that the Reynolds number and the Peclet number 
are sufficiently large so that mixing is dominantly due to turbulence. 

Integral equations are obtained by integrating each of equations (EH) over a cross- 
section of the plume. We define a surface r = b{z,t) representing the boundary of the 
plume over which entrainment of ambient fluid into the plume occurs and impose the 
boundary condition on this surface. 


u{b,z,t) 


dt 


+ w(6,z,<)— -h 


( 2 . 2 ) 


where the entrainment velocity, Ue, is the radial velocity of the ambient fluid across 
r = 6(z, t) (note Ue < 0 when ambient fluid is entrained into the pl u me). W e note 
that this approach differs from that taken by Craske fc van Reeuwiik ( 2015ai r^ who 
adopt the characteristic length and velocity scales defined through the cross-sectionally 
averaged fluxes of mass and axial momentum which themselves are defined in terms of an 
additional length scale at which the mean axial velocity can be considered negligible in 
contrast to the mean axial velocity at the centreline. As turbulent plumes typically exhibit 
non-uniform radial profiles of axial velocity and concentration of species, with decaying 
velocity and concentration as the radial distance from the plume axis, in laboratory 
experiments or direct numerical simulations there is a choice as to how to define the 
plume edge, taking, for examples, the radial distance at which the axial velocity or 
concentration reach a specified threshold, or the radial distance at which the mass flux 
included in a cross-sectional integral captures a specified proportion of the total mass 
flux induced by the plume. We discuss below how these choices are represented in our 
model. 


Integration of the point-wise conservation equations EH over a plu me cross-section , 
using Leibniz’s theorem for interchanging differentiation and integration (jFlanderslll97,‘l|) 
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together with the boundary condition (12.21) . gives 

Jb 9 
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r dr 
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/ rg!^w'dr-bg'^{b,z,t)ue-bu'g'^{b,z,t)-bw'g!^{b,z,t)-^, (2.3c) 


9z 


9z’ 


representing conservation of mass, momentum and buoyancy, respectively. 

We define the integral volume flux as Q = 6^W (note, under the Boussinesq assump¬ 
tion, Q can also be referred to as the (specific) mass flux), the (specific) momentum flux 
as M = b'^W^ and the buoyancy flux as F = b^Wg'. Here W and G' denote the cross- 
sectionally averaged mean axial velocity and reduced gravity of the plume, respectively, 
and are given by 


W = 


2 

If 


rw dr, and G' = 


2 

If 


rgr dr. 


(2.4a, 6) 


We then have 


pb pb 

/ ruJ^ dr = iS'M, and / rw^dr = ^(j)F, (2.5a 

Jo Jo 

respectively, defined as 

f r {w — W)"^ dr, 

Jo 


b) 


10 Jo 

where S and cj) are momentum and buoyancy ^shape factors' 

f.b 

5 = 1 + 


621+2 


and 


= 1 + 


- j\{w-W) {Tr-C) 


dr. 


( 2 . 6 ) 


(2.7) 


IFWg' jQ 

Thus, the shape factors quantify the effect of non-uniform radial profiles of the axial 
velocity and density on the rates of transport of momentum and scalar species. Crucially, 
we note that 5^1, which is representative of the more rapid transport of momentum 
than mass in the plume, unless the mean axial velocity is radially invariant within the 
plume (in which case 5=1 and ^ = 1). 

The plume edge, given by the surface r = b{z,t), is chosen to be at a radial distance 
such that the boundary terms in (12.31) are negligible in comparison to the integral terms 
(e.g. W(b) WjU'w'{b) <C W'^, etc.). Furthermore, the entrainment assumption of 
Morton et all 1 1956h allows the entrainment velocity to be written as Ug = —kW where 


k is the entrainment coefficient. This simplification of the turbulence mixing dynamics 
assumes a similarity of turbulent structures at each height in the plume. The entrainment 
coefficient must be determined empirically and, for fully developed plumes far from the 
source, a con stant value k = 0.1 + 0.01 is appropriate ( Morton et al\ 1956t List 1982 : 
Woodd 2 OIOII . However, near to the source (or for non-ideal initial conditions), the en¬ 


trainment coefficient may vary substanti ally as the flow develops (jWang & Lawl 12002 ; 
Kaminski et a,l. 20051 : Ezzamel et aLlEoibh . In this study we take a constant entrainment 
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coefficient (with the exception of appendix |A] where we also examine the steady solutions 
of a model with an entrainment coefficient that varies as the flow develops). 

We note from (1^ that the value of the momentum shape factor is tied to the choice of 
the plume width, b[z,t). If we assume the mean axial velocity has a self-similar Gaussian 
profile in the radial direction, as inferred from labor atory experiments on steady plumes 
I Papanicolaou fc Lisdll988 : Shabbir fc Georg jll994ll . so that the mean axial velocity of 
the plume can be written as 


w (r, z, t) = W{z, t)e 


-A jB? 


( 2 . 8 ) 


where W(z, t) is the mean axial velocity on the plume axis and R{z, t) is a characteristic 
length scale for radial variation, then the shape factor is given by 


S = 


52 ('i + 

2i?2 


(2.9) 


The plume width can be related to the length scale of the Gau ssian radial profile R {z, t) 
through a threshold on the mean axial velocity. For example, Morton et al. ( 1956l l and 
Papanicolaou fc List ( 1988h define the characteristic length scale of a Gaussian plume as 
the radial position at which the axial velocity is a factor of 1/e of the centreline value, 
and this results in 5 = 1.08. 

The entrainment hypothesis and use of shape factors for the integral fluxes of momen¬ 
tum and buoyancy allow the integral conservation equations (12.3p to be written in terms 
of the integral fluxes of volume, momentum and buoyancy as, 

91 ' 

M 


0 

dt 


+ ^ = 2fcyM, 

az 


(2.10a) 




0 

dt 


QF 

M 


dz po 


2rw'^ dr, 


dt dz 


( 2 . 10 &) 


I; 


2rg'^w'dr. (2.10c) 


While the system (12.101) applies to a spatially and temporally varying ambient den¬ 
sity field, for the remainder of this study we consider plumes in an unstratified ambi¬ 
ent, with Pa constant. The shape factors S and (j) could be spatially and temporally 
varying, as the radial profiles of the axial plume ve lo city and reduced gra vity evolve 
( Carazzo et al\ 200d Craske fc van Reeuwiik 2015al f^: lEzzamel et al\ 2015ll . However, 
for fully developed plumes, laboratory experim ents suggest that the axial plume velocity 
and reduced gravity attain se l f-similar forms dMorton et al. 19561 : Papanicolaou fc List 
1988t Shabbir fc Georgjll994 lEzzamel et all 2015 ) such that the shape factors can be 
taken to be constants. Here we consider the simplest case of constant S and (j)= 1. 

In the integral equations for conservation of momentum (I2.10b|) and conservation of 
buoyancy (I2.10cl) we have retained integral terms that represent turbulent axial diffusion 
of momentum and buoyancy, respectively. Typically, in steady plume models these diffu¬ 
sive terms have been neglected (we refer to the system of integral eq uations (12.101) w ithout 
the turbulent diffusive terms as the ‘non-diffusive’ system). Indeed. iMortonI ( 197ll) argues 
that 

3 1 3 3 13 

— (ru'w') , and — (w'gA) — (ru'gA) , (2.11) 

dz \ / r dr dz r dr 

and therefore the contributions of the turbulent diffusion terms to the plume dynamics 

are of a similar magnitude to boundary terms that are neglected. We note some studies on 
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momentum driven jets from maintained sources retain the axial der i vatives of quadratic 
fluctu ation terms (see e.g. Shabbir fc George 1994 IWang fc Law! I 2 OO 2 I : lYannapoulos 
2006li . When these turbulent diffusion terms are neglected, the steady integral equa¬ 
tions with pure plume boundary conditions (Q = 0 , M = 0, F = Fq at z = 0) have 
well-known power-law solutions ( Morton et al. 1956ll . and these solutions are little al¬ 
tered by the momentum shape factor; for S' ^ 1 the steady solutions of the non-diffusive 
system of equations are 


Q = Qo{z) = qoz 


— n r5/3 


M = Mq = 


F = Fn 


( 2 . 12 ) 


qo = 


5 Vio/ 


(2.13) 


where 

1/3 

In the limit S —^ 1 the solution of iMorton et al\ (Il956ll is recovered. Furthermore, the 
effective radius of the plume, bo{z) = Qq/^/Mq — Qkzjb, is independent of the shape fac¬ 
tor. Therefore, the momentum shape factor cannot be determined from measurement of 
the plume radius alone, in contrast to the entrainment coefficient that can be determined 
u sing the spreading rate of the steady plume. 

Scase fc Hewitt ( 2012[ l advocate modelling of the turbulent diffusion terms (partic¬ 


ularly the diffusion of momentum) in (12.101) in the unsteady plume mo del to obtain a 
well-posed system of equations. Recently, Craske fc van Reeuwiik (201530) have demon¬ 
strated the relatively weak role of turbulent diffusion for momentum driven jets, but that 
diffusive effects in the jet occur due to the dep arture of flow variables from se lf-similar 
forms (type II dispersion in the nomenclature of Craske fc van Reeuwiikll2015al f3l . How¬ 
ever, we suggest that the dominant dynamics for turbulent plumes can be described by 
the non-diffusive system of equations. Indeed, we show below that the inclusion of dif¬ 
fusive terms modelling turbulent diffusion can lead to difficulties in the integral model. 
We therefore analyse the non-diffusive system and show that a momentum shape factor 
that differs from unity is sufficient to obtain a well-posed model. 


3. Difficulties associated with the turbulent diffusive terms 

In an analysis of a non-diffusive unsteady plume model with shape factors S = 1 and 
(j) = 1, corresponding to top-hat radial profiles for the axial velocity and reduced gravity 
(i.e. w and ^ are radially invariant for r ^ b and equal to zero for r > b), IScase fc Hewitt 
( 2ni2h assess the well-posedness of the system of equations by introducing small harmonic 
perturbations to the source buoyancy flux and examine the growth of the perturbations 
downstream of the source. A linear analysis shows that the perturbations grow with dis¬ 
tance from the source, indicating instability, and, more importantly, the growth rate of 
the perturbations increases without bo und as the frequency of the harmonic oscillation 
of the source buoyancy flux increases ( Scase fc Hewittl[20I2ll . The non-diffusive system 
of equations with S' = 1 and (j) = 1 are therefore ill-posed, as t here is a loss of con¬ 
tinuous dependence of the solution on the boundary conditions ( Joseph fc Saut 1990l : 
Scase fc Hewitdl20I2ll . The ill-posedness is manifest in numerical solutions of the system 
of equations by an inability to compute solutions that are independent of the truncation 
implicit in the numerical scheme (e.g. grid-scale dependence for finite-difference meth¬ 
ods). 

The ill-posedness identified bv I Scase fc HewitO ( 2012 ) in the non-diffusive system when 
top-hat profiles are assumed (i.e. when S' = 1) may be due to a missing physical process 
that provides a mechanism to curtail the unbounded growth of the arbitrarily short- 
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wavelength/high frequency modes iTosenh fc SautlllOOC : 


tempt to regularize the ill-posed system, IScase fc Hewitd (j2012l l introduce a phenomeno- 


Scase fc Hewitt 2012[l . In an at- 


logical model for the diffusive term in the momentum balance (I2.10&|) . appealing to 
Prandtl’s mixing length theory for turbulent eddy diffusion, 


0 

Vz 


f’’ 

K r, 0 

01T 

/ 2rw' dr r 

— b^— 

bW - 

Jo 

2k dz 

dz 


(3.1) 


where k > 0 is a dimensionless parameter characterizing the diffusion of momentum 
through the action of turbulent eddies whose length scale is set by the rad ius of the plume. 
The e quation expressing the balance of axial momentum proposed by IScase &: Hewitt 
(I 2 OI 2 II is then given by 


(3.2) 


dt dz M 2k M dz dz \Q 

Scase fc Hewitd ( 2012h provide numerical evi dence that the dif f usive term leads to a well- 
posed system of equations. Note that, while IScase &: Hewitd ( 2012 1 take S' = 1, in the 
a nalysis presented below we analyse the more general problem of S ^ 1. 

Scase fc Hewitt ( 2012ll show that steady power-law solutions for pure plume boundary 


conditions exist for the diffusive system (with k > 0) and the str uctural change to the 
system of equations leads to a small modification of the solutions of Morton et al ] (fl^ 
(given bv 12.121 with S = 1). The steady solutions of t he ‘regularized’ system of equations 
with turbulent diffusion of momentum are given by (IScase fc Hewittll2012r) 


= 90 (1 - 


IQS'/ 


-1/3 


,5/3 






(SH) 


= Fn. 


(3.3o) 

(3.36) 

(3.3c) 


However, we find that these power-law solutions are spatially unstable. Indeed, taking 
perturbations of the form 

(z) = ( 2 ) (1 + egf(^)) , ( 2 ) = (^) (1 + eMf (z)) , 

(3.4a, h) 

(noting that the buoyancy flux F{z) = Fq for a plume in an unstratified ambient) where 
e > 0 is an ordering parameter, and linearizing the governing ordinary differential equa¬ 
tions (for e <C 1) we obtain. 


- jM.®"' = 0 . (3-5») 

55'^ dz2 dz2 j ^ IQS'/ ^ dz 

+ 5<?r+4(^-5)Mf"> = 0. (3.56) 


Seeking solutions of the form g) = giz“ and = miz“, we find non-trivial 

solutions are possible ifa = Q!i = —I or 

55' -b 4 k ± V25S'2 -h 240 S'k - 4^2 

a = a± = ---. (3.6) 

OK 
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K K 

Figure 1. Spa t ial ar owth rates Re (a) as a function of the diffusion coefficient for (a) the 
IScase Sz HewittI ll2012f) system with momentum diffusion (13.211 . (b) the aiternative form of mo¬ 
mentum diffusion (I3.8II . and (c) diffusion of buoyancy (13.1111 . In (a) and (b) different values of 
the momentum shape factor are shown, with 5 = 1 (solid fine), S = 1.5 (dashed fine) and S = 2 
(dash-dot line). In (c) the three branches of the growth rates are distinguished, with one branch 
corresponding to a real mode (solid line), and two branches that form a complex conjugate pair 
(dashed and dash-dot lines). 


At least one of a-i- > 0 for k > 0 for any S 1 (figure [T^) and therefore the steady 
solntions of the IScase fc HewittI (l2012i) ’regularized’ model are spatially unstable. We 
note that, in the limit k —>■ 0, we find a+ no longer appears in the analysis while a_ —>■ 
— 10/3 and therefore the steady solutions be come spatially stable (as expected, since the 
‘regularized’ system reduces to the classical [Morton et al. 1 19561) plume model). Spatial 
instability of the steady solutions is also found in a stratified ambient (not shown here) 
and in this case leads to severe difficulties as steady solutions cannot, in general, be 
found analytically and the spatial instability precludes the nu merical computa t ion o f 
steady solutions. Therefore, while the diffusive term introduced bv IScase fc HewittI ( 2012 ) 
into the momentum balance resolves the ill-posedness in the unsteady plume model, the 
modification of the governing equations renders the steady solutions spatially unstable. 
While our analysis here considers pure plume boundary conditions (i.e. Q(0) = 0, M(0) = 
0 and A(0) = Fg) in appendix [B] we show that spatially stable steady solutions with 
arbitrary initial conditions are not possible. 

An alternative (although similar) form of the diffusion term can be formed by re¬ 
ordering the cross-sectional integration following the mixing-length parameterisation of 
the fluctuation vertical momentum flux (with boundary terms introduced that are as¬ 
sumed small and subsequently neglected), so that we write 




b'^W- 


dW_ 

dz 


(3.7) 


in place of dSU). We then obtain the following expression for conservation of momentum, 

(3.8) 


dQ 9M _ QF K d 
dt dz M 2k dz 


r 9 

fM\] 

Vm dz 



The corresponding steady solutions are 

-)(SH2) _ 


-1/3 


Q/-^i=gg(l + A)-“,5/3^ (3.9a) 


= mg 




-2/3 


A/3 


F, 


(SH2) 


= Fn 


(3.96) 

(3.9c) 
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and a spatial stability analysis gives growth rates a = ai = — 1 or 


a = a± = 


55 - 3k ± V2552 + 1705k + 49^2 
Iok ■ 


(3.10) 


Thus, as with the IScase fc Hewittl ( 2012 1 momentum diffusion term, there is a growing 
mode (figure [T}d) and the steady solutions are spatially unstable. 

Finally, we consider the change to the classical steady solutions of Morton et al ] (Il956h 
that occurs when a diffusive term is included in the equation for conservation of buoyancy 
(|2.10cl) rather than in the equation for conservation of momentum (i.e. we now take 
K = 0). We again appeal to the Prandtl mixing length theory to model the fluctuation 
of the flux of buoyancy and obtain the following diffusive equation for the conservation 
of buoyancy. 


0 

dt 


QF\ 
M ) 


dF 

dz 


Ki 0 
2k dz 




(3.11) 


_y/M dz 

where ki is a dimensionless parameter that describes the strength of the diffusion of 
buoyancy. The classical steady solutions (12.121) for pure plume boundary conditions re¬ 
main as solutions of the system (and other possible solutions cannot satisfy the pure 
plume boundary condition at 2 = 0). However, a spatial stability analysis of the steady 
solutions (here including perturbations to the buoyancy flux in addition to the mass and 
momentum fluxes) shows that small amplitude perturbations to the steady solutions are 
amplified unless the magnitude of the turbulence induced is relatively large (ki > 4/3, 
figure [T}:). 

The consequences of these analyses are that, while axial diffusivity affects may capture 
some important unsteady processes in the mixing of momentum or buoyancy, the form 
of the parameterization significantly alters the stead y states attainable b y the system. 
Indeed, the states equivalent to those established by Morton et ol\ ( 19561 ) have become 
spatially unstable and therefore the unsteady model with axial diffusivity is not able to 
describe the steady states. 


4. Well-posedness of the hyperbolic unsteady plume model 

As the turbulent diffusive terms modelled phenomenologically using a Prandtl mixing- 
length approach leads to a pathology in the in tegral plume model whereby the well- 
established steady states of iMorton et al\ 1 1956ll cannot be obtained, the regularisation 
through eddy diffusion, while potentially physically appealing, is mathematically prob¬ 
lematic. We therefore seek an alternative regularization of the system of unsteady integral 
equations; specifically, we examine the non-diffusive s ystem of equat ions (i.e. we neglect 
the turbulent diffusive terms in 12.101 as suggested by iMorton 197lll but allow the mo¬ 
mentum shape factor to differ from unity to describe the different rates of transport of 
mass and momentum that results from non-uniform radial profiles. Thus, from here on 
in we study the non-diffusive system of equations. 




dQ ^dM _ QF 0 dF 


dt 


dz 


M 


dt \ M J dz 


-I- — = 0. (4.1a, b, c) 


The system 63 ) is strictly hyperbolic when 5 > 1, with three real characteristics wave 
speeds given by 

Co = M/Q, and c± = ^5 ± v'5'(5- 1)) M/Q, (4.2) 

We note, however, in the limit 5—^1 there is a loss of hyperbolicity since the eigenvectors 
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of the characteristic equation for the system dUD do not span althou gh the wave 


speeds remain real and equal to co; the system is formally parabolic for 5" = 1 (jScase et al 


I 2 OO 66 II . As we demonstrate below, the change in the characteristic structure, and so 
change in type, of the governing equations that occurs for 5" = 1 fundamentally alters 
solutions of the system of equations. (Note, if a shape factor for the buoyancy flux that 
differs from unity is included in 14.1b . then cq = 4>M/Q while c± remain unchanged, and 
a loss of hyperbolicity occurs if </> = ^ ± S {S — 1).) 

In order to establish well-posedness of the system of equations, we analyse the evolution 
of small perturbations to the steady solutions. It is convenient for subsequent analysis and 
numerical computations to factor-out the steady solution and to introduce a new spatial 
coordinate, x = {qQ/mQ)z^^^. We consider the stability of the steady solution to small 
perturbations. Therefore, we introduce an ordering parameter e <C 1 and perturbations 
to the steady solution of the form. 


Q{x, t) = Qo{x) (1 -b eq{x, t )), 
M(x, t) = Mq(x) (1 -b em{x, t)) 
F{x,t) = Fo {l + ef{x,t )), 

and linearize the governing equations dUD to obtain 


9t 


3 dx 


C-q 

X 


(4.3a) 

(4.36) 

(4.3c) 


(4.4) 


where q = {q,rn, f )^ and the matrices A, B, and C are given by 


A = 



B = 



and 



When S' = 1 the system (14.11) was shown to be ill-posed by Scase fc Hewitt! ( 20121) who 
examined the response of the linearized system of equations to small amplitude harmonic 
variation of frequency w in the source buoyancy flux, and found a closed form solution 
for the perturbations in terms of special functions. It was shown that the amplitude of 
the perturb ations grows with dis tance from the source as exp {vTSwx/d} 

for a; » 1 ( Scase fc HewittI 20121) . Thus, the steady solutions are linearly unstable, as 
the amplitude of the perturbations grows as they propagate away from the source and, 
furthermore, the system of equations (14.11) with S' = 1 is ill-po sed since high frequen cy 
fluctuations increase in amplitude most rapidly with no cut-off (|Scase fc Hewittll2012l) . 

For S > 1 no closed form solution of the linear system (lOl) can be obtained. We 
therefore examine the evolution of perturbations in the far-held through an asymptotic 
analysis of the linear system (14.41) in the neighbourhood of the irregular singular point 
X —>■ 00 ((Bender fc Orszad 119781) . We consider two stability problems; an initial value 


problem where the evolution of an arbitrary initial perturbation with compact support. 


q(x, 0) = qgix) with qgix) —>■ 0 as x —>■ 00 , (4.6) 

is investigated, and a boundary value problem where a fluctuation at the source x = 0 
is imposed and the response of the system downstream of the source is examined. Both 
of these problems can be analysed conveniently through the use of integral transforms 
of the linear system (14.41) : a Laplace transform in time for the initial value problem 
and a Fourier transform in time for the boundary value problem. Denoting the Laplace 
transform of q{x,t) as q{x,p) (with the Fourier transform obtained by taking p = iuj 
where uj is the frequency of the harmonic source fluctuation imposed in the boundary 
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value problem) we obtain, 

.. 4 dq „ 

pAq + 3®-^ + 


(4.7) 


The f ar-field behaviour is obtained conveniently by letting q = ( Bender fc Orsza^ 

1978h with g{x) a singular function and f(x) regular as a: —>■ oo, so f(x) = /q + + 


f^x ^ + ... for X » 1. The linear system (IT71) can then be written as 


4od/ 


3 dx 


4 o dg , 1 


Ti ® — f 1 ~b ~ 6 —— I — C 1 f — 0, 


3 dx X 


(4.8) 


A leading-order balance requires g{x) ^ pXx as x —>■ oo. Then A and /g satisfy the 
generalized eigenproblem A/q = —4A6/g/3 and therefore, 

3 

1 


A = Ao = --, with /(, = /oo = (0,0,1)^, /^ = /^g = (-1,0,1), 


(4.9a) 


A = A+ =-- 


3 / ^+^7*5 (5-1) 


5 


with /g = /g = 1, 


VS{S-l) 


,1 


= /^+= (5+v'5(5-l),-l,0), 




1 - 


with /g = /g_= 1 ,-- 


(5+v^5(5-l)) ^ 

VS{S-1) 


(4.96) 

T 


= ft = {S- v'5(5-l), -1, o) 


(4.9c) 


where /g denotes the left eigenvector satisfying f^A = —4A/o 6/3. 

To proceed further we let g{x) ^ pXx + fi log x for x » 1. Substitution into (14.81) and 
balancing coefficients of x gives, at order 0(l/x), 


7'('4-I--A6)/]^-)-( —p,B -f C I /g — 0 


(4.10) 


and so, by multiplying on the left by ft we hnd 


3/^C/g 
^ 4/^e/g’ 

(4.11) 

and therefore 


/r = p,o = 0 when A = Ao, 

(4.12o) 

13 5(25-1) ^ , 

a = = -, = when A = A+, 

8 16v'5(5-1) 

(4.126) 

13 5(25- 1) , , , 

u = u =-1-, = when A = Aj.. 

8 16v'5(5-1) 

(4.12c) 


The coefficient vectors (/j^, / 2 , etc.) can be determined from higher order terms in the 
expansion but are not given here. 

The leading order behaviour in the far-held is therefore given by, 

q{xA) ~ foX^5 {t - Ax), 


(4.13) 
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Figure 2. The exponents of the far-held algebraic growth rates of the amplitude of the volume 
hux perturbation, p.+ (dashed line) and fi- (solid line), as functions of the shape factor S. The 
perturbation to the volume hux grows with increasing distance from the source when /r_ > 0 
which occurs for S < 25/24. 


corresponding to the propagation of discontinuities whose strength varies algebraically 
with distance from the source. The amplitude of the perturbations grows if 1 < S' < 25/24 
whereas the amplitude decays algebraically if S > 25/24 (figure [2]). Importantly, the 
algebraic growth rates fa do not depend on the transform variable p. Therefore, for the 
initial value problem the evolution of perturbations in the far-held is independent of the 
spatial length scale of the initial perturbation, whereas when considering a boundary 
value problem the growth rate of the perturbations do not depend on the frequency of 
the harmonic boundary oscillations. Therefore, the far-held asymptotic analysis shows 
that the system 63 ) with S > 1 is well-posed as there is a continuous dependence of 
the solutions on the initial or boundary data and small scales (either spatial scales in 
the case of an initial value problem or temporal scales for a boundary value problem) 
are not amplihed more rapidly than longer scales. This is in contrast to the far-held 
limit o f the solution of a boundary value problem that is determined by Scase fc HewittI 
( 2012h . where exponential growth of the amplitude of perturbations are found with a 
growth rate that increases exponentially with -y/w, and therefore when S' = 1 the system 
of equations is ill-posed. Numerical solutions demonstrating the evolution of a localised 
initial perturbation are presented in (|5] below. 


5. Solutions of the well-posed unsteady plume model 

We consider now solutions of the nonlinear system of equations 63 with S > 1. As the 
system of equations 63 is hyperbolic in this parameter regime, solutions may exhibit 
discontinuities, across which we enforce the following jump conditions that respectively 
impose the conservation of mass, momentum and buoyancy huxes over a discontinuity 
at z = Zs{t), moving with velocity c = dzs/dt, 


Q — c 


91 

M 


= 0 , 


SM — cQ 


= 0 and 


F-c 


QF 

M 


= 0. (5.1) 
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These jump conditions admit two different types of discontinuous solutions, which are 
more easily analysed in terms of the primitive variables b, W and G'. Further denoting 
the values of the dependent variables either side of the discontinuity with superscripts 
and corresponding to z = and z = z~, respectively, we find that provided the 
velocity is discontinuous (VF+ ^ W~) then the reduced gravity G' is continuous (i.e 
G"'*' = G'~). Furthermore in this case, by eliminating 6+ and b~, we find that 

(W+ _ W-) (c^ - cS{W+ + W-) + SW+W-) = 0. (5.2) 

When the shape factor is equal to unity, the only solution is c = VF+ = W~ and this 
contradicts the assumption that there is a discontinuity. Thus for S' = 1 it is not possible 
to construct discontinuous solutions. However, when S > 1, we find 

c = \S{W+ + W-) - i (S2(VF+ + W-f - 4.SW+W-f'^ , (5.3) 


where the negative root has been chosen so that W~ > c, a condition required for 
causality. The other type of discontinuity occurs when W~^ = W~ = c. Then the radius 
of the plume is also continuous, b'^ = b~, but the reduced gravity may be discontinuous 
and potentially even unbounded. We refer to this latter form as a ‘contact’ discontinu ity, 
which occurs due to a linearly degenerate field in the system of equations ( Laxlll97il . 

To calculate numerically solutions of the n onlinear hyperbolic sys t em of equations, we 
employ the non-oscillatory central scheme of Kurganoy fc Tadmor ( 2000ll with a third- 
order total yariation diminishing Runge-Kutta time stepping scheme ( Gottlieb fc ShiJ 
Il998ll using an adaptiye time step that ensures the CFL number remains fixed at 1/8 
to maintain numerical stability. The high -resolution central scheme i s a shock-capturing 
numerical method for conseryation laws ( Kurganoy &: Tadmor 2000ll and has been used 
extensiyely to compute solutions of nonlinear hyperbolic systems. For numerical conye- 
nience, we factor-out the steady solutions from the dependent yariables and compute 
solutions in the transformed spatial coordinate x as, although the system is not au¬ 
tonomous in this representation, the characteristic waye speeds remain bounded near to 
the source (i.e. as x —> 0) when the steady solution is realized, whereas the waye speeds of 
the system diyerge as z —> 0 and therefore small time-steps are required to maintain 
numerical stability. 

We consider first the eyolution of a perturbation to the steady solution (12.1211 to 
confirm the far-field asymptotic analysis in 0 We take as an initial condition a Gaussian 
perturbation (in the transformed coordinate x) to the steady momentum flux of the form 


M(x, 0)/Mo(x) = 1 -f 0.01e-^°(‘"-®)' 


(5.4) 


while the yolume flux and buoyancy flux are not perturbed from the steady values, and 
take the shape fact or to be S = 1.1. We compu te numerically the evolution of the per¬ 
turbation using the Kiirganov fc Tadmor ( 20nol l central scheme for the nonlinear system 
rather than using the linearised equations since the linearisation introduces a subtle 
structural change in the governing equations; for the linearised system, all field are ( triv¬ 
ially) linearly degenerate and so discontinuities are of the contact discontinuity type (jLax 
197,3[) whereas two fields i n the nonlinear system are genuinely nonlinear. We find that 


the Kurganov fc Tadmoil ( 200(TI 1 scheme resolves accurately the shocks and rarefactions 
associated with the genuinely nonlinear fields, but has l ess accuracy when contact dis¬ 
continuities occur for linearly degenerate fields (see e.g. iKurganov fc Petroval 1200011 . In 
order to track the evolution of the perturbations to long times (and so large distance 
from the source), we implement a moving numerical domain, using the characteristic 
wave speeds of the linearized system to advect the lower and upper grid points. The 
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numerical domain is then advected downstream and the number of grid points increases 
as the spatial extent of the perturbation grows in order to maintain a specified numerical 
resolution. 

In figure[3]we illustrate the evolution of the initial perturbation for dimensionless times 
t ^ 40, showing the perturbation to the plume radius, mass flux and buoyancy flux as 
functions of the scaled spatial coordinate x. The initial perturbation develops into a pulse 
whose spatial extent grows as it propagates. The amplitude of the perturbation slowly 
decreases, demonstrating linear stability of the system of equations for S' = 1.1. The three 
discontinuities are apparent in the evolving perturbation, with the contact discontinuity 
that propagates at the speed of the intermediate characteristic most clearly seen in the 
perturbation to the reduced gravity, G'. 

We introduce integral measures of the amplitude of the perturbations, with 



(5.5) 


for the volume flux, and similarly for the fluxes of momentum and buoyancy. Figure 
Hi shows the time evolution of I{Q), I{M) and I{F) when the shape factor S = 1.1. 
Following a short transient reorganization of the initial condition for t < 5 (not shown 
on the figure) the perturbations to the steady solution decay as they propagate through 
the domain. The perturbations to the steady fluxes of volume and momentum closely 
follow the prediction of the far-held analysis for t > 100, at which time the signal of 
the perturbation that travels with the fastest characteristic has reached x ~ 197. The 
decay in the perturbation to the buoyancy flux is less rapid and the rate of decay is 
diminishing as time progresses. This is expected from the far-held asymptotic analysis of 
the linearised system where a component of the buoyancy hux is found to be advected 
at the speed of the plume without change in amplitude (see 14. 12 all . 

For S < 25/24 our far-held asymptotic analysis shows that the steady solutions are 
linearly unstable. Numerical solutions of the governing equations support the asymptotic 
analysis, as illustrated in hgureHr which shows the time evolution of the integral measures 
of the size of the perturbations for a shape factor S = 1.02. For early times {t < 10) the 
perturbations to the steady solutions grow algebraically with a growth rate that is close 
to the growth rate predicted from the far-held analysis. We note that at early times the 
perturbation has not propagated far downstream from the localised initial condition and 
therefore the far-held asymptotic analysis would not be expected to describe fully the 
evolution. At later times {t > 30) the growth of the perturbations deviate substantially 
from the predictions of the linear far-held analysis as nonlinear effects begin to inhuence 
the evolution. 

Numerical solutions of the governing equations also support the far-held asymptotic 
analysis of the linearized system when a harmonic oscillation of the source boundary 
conditions are imposed. Figure [5] illustrates the spatial structure of the perturbation to 
the steady volume hux at time t = 99 for a shape factor S = 1.05 (hgure[5K) and S = 1.03 
(hgure[S]D). The far-held analysis predicts linear stability when S = 1.05 and instability 
when S = 1.03, and this is observed in the numerical solutions with the amplitude of 
perturbations decaying with increasing distance downstream when S = 1.05 (hgure[5^), 
in contrast to the growth in the amplitude of the perturbations that is seen when S = 1.03 
(hgureEb). 
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{Q/Qo - 1) X 10® (M/Mo - 1) X 10® {F/Fo - 1) x 10® 



Figure 3. The perturbations to the fluxes of volume, QjQo — 1, momentum, M/Mo — 1, and 
buoyancy, F/Fo — 1, and to the plume radius, b/bo — 1, the plume velocity, W/Wq — 1, and 
the reduced gravity, G' IGq — 1, as functions of the scaled distance from the source, x, for 
dimensionless times t = 0 (orange, colour online), t = 10 (red, colour online), t = 20 (green, 
colour online), t = 30 (blue, colour online), and t = 40 (black, colour online). The momentum 
shape factor is fixed at a value 5 = 1.1. The numerical solution scheme adopts a moving and 
expanding spatial grid with a fixed grid spacing of Aa; = 0.01. 
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Figure 4. Integral measures of the size of perturbations from the steady solutions for the 
volume flux, T{Q) (solid line), momentum flux, I{M) (dashed line), and buoyancy flux, T{F) 
(dotted line) as functions of the dimensionless time t, for a shape factor (a) S = 1.1 and (b) 
S = 1.02. The thick grey line illustrates the predicted growth rate of perturbations of 
with (a) ~ —0.49 for S = 1.1 and (b) fj,+ ~ 0.65 for S = 1.02 obtained from the far-held 

asymptotic analysis (14.12611 . 


6. Similarity solutions for the unsteady evolution of plumes following 
an abrnpt change in the source buoyancy flux 

We examine now the nonlinear evolution of plumes following an abrupt change in the 
source conditions. Specifically, we consider solutions of the system of nonlinear evolution 
equations with S > 25124: (so that the steady states are linearly stable) for an initial 
boundary value problem in which the magnitude of the buoyancy flux at the source is 
instantaneously changed from Fq to at t = 0. For t < 0, the plume is in a steady state 
given by ()2.12p and ()2.13l) . We calculate the unsteady evolution as the plume adjusts to 
the new steady state in which the buoyancy flux, Fq, in ()2.12l) and (12.131) is replaced by 
Fi. Thus, the initial conditions are given by (I2.12|) and (12.131) . while for t > 0 the new 
source conditions for the plume are given by F{0,t) = Fi and Q{0,t) = M{0,t) = 0. 

The numerical solutions demonstrate that the adjustment occurs by the advection 
of an unsteady ‘pulse’ through the environment (see figures [SI [5] and (TU] for examples 
of computations), which is modelled by the nonlinear evolution equations (14.11) . In this 
section we draw out the self-similar adjustment that occurs in the dependent variables. 

The similarity variable is established through the following scaling analyses, which 
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-0.04 0 0.04 

QIQo — 1 

Figure 5. The perturbation to the steady volume flux due to harmonic oscillation to the source 
buoyancy flux with period 27r as a function of the scaled spatial coordinate x at dimensionless 
time t = 99 for a shape factor (a) S = 1.05 and (b) S = 1.03 . The red dashed line (colour 
online) illustrates an algebraic growth of perturbations of the form cx^ where (a) g = —0.12 for 
S = 1.05 and (b) g — 0.26 for S = 1.03. The prefactor c is selected to provide an envelope of 
the numerical solution. 


require all of the terms in governing partial differential equations to be of the 

same order. To this end, to balance all the derivative terms we require M ^ Qz/t. 
Turning then to the ‘source’ terms, from (14.Ih ). we have Qjz while from (14.1b ) 

Mjz ^ QFi/M, where we have used F ^ Fi a.s the scale of the buoyancy flux. Together 
these yield z* ^ Fit^ and we therefore seek similarity solutions to the system in terms 
of this gearing between the spatial and te mporal scales . We no te that t he existence 
of th is similarity grouping was identified by Scase et ^ ( 20066h (see also Scase et al 


20081) . although they did not construct the complete similarity solution; indeed for their 


system with 5” = 1 it was not possible to do so because the ill-posedness of the system 


manifes ts itself as irreconcilable divergences in the solution (see Q . Instead IScase et al 


(j 2006 &l l found a separable solution which did not satisfy the boundary conditions, but 
did capture some of the features found in their numerical computation. We discuss in 
appendix |D] the counterpart to their separable solution in the now regularised system of 
governing equations. Before analysing the form of the similarity solutions, we note that 
there is another reason for anticipating the similarity scaling deduced above: the system 
is hyperbolic for 5 > 1 and all the characteristics are proportional to W = M/Q oc 
The characteristic curves are of the form dzj dt oc W, thus also admitting 
the similarity scaling z^ ^ Fit^ and it is through the evolution of the characteristics that 
the system adjusts to its new state. 

We now seek similarity solutions for the fluxes of volume, momentum and buoyancy, 
which we write in the following form 
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where the similarity variable is given by 


and F = FiF{ri), 

( 6 . 1 ) 
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and the similarity functions Q(? 7 ), M{vi) and F{rf) are to be determined. In this form 
the steady state given by (12.121) corresponds to constant values of Q, M and F. It is 
convenient to further substitute Q = rjq^ M = rf'rh and F = rf f because this simplifies 
the governing equations, which are now given by 
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Symbolically, we write this coupled differential system as Dr]dq/dr] = 6, noting that 
both the matrix D and the vector b are only functions of the depe ndent variables q = 
{q, fh, /). We note that the ‘separable’ similarity solutions derived bv IScase et all ( 20066 1 
correspond to q = {q,fh,f) = constant and these are discussed in appendix [D] For the 
motion driven by an abrupt change in the magnitude of the buoyancy flux at source, the 
vital parameter is the ratio of the initial to final buoyancy fluxes at the source given by 
F = Fq/Fi . We note that the initial conditions correspond to 

l'jrl/3 j-2/3 


(q, m, f) = 




V 




(6.4) 


In terms of the similarity functions, these are enforced in the far-held {q —> oo) and 
correspond to the region that is unaffected by the change of the buoyancy hux at the 
source. The source conditions demand that 


iq,rnj) 


1 1 1 
q^ ^2’ jj3 


as 
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(6.5) 


Constructing the similarity solutions then corresponds to integrating the governing sys¬ 
tem of ordinary differential equations (16.31) . subject to these conditions. 

The matrix D is singular when 


m 


16 

y 


3 fh 


= 0 . 


( 6 . 6 ) 


This corresponds to locations where q/fh = 4/3 and q/fh = 4 (S' ± \/S^ — S) /3. These 
are singular points of the governing system of equations and are of signihcance because 
the dependent variables may have discontinuous gradients at these locations. We show in 
appendix [C] how to derive the local behaviour close to the singular points; this is required 
to initiate numerical integration from these locations. 

The similarity solutions may also feature discontinuous solutions, in which case we 
define the shock position, Zs{t), scaled according to the similarity variables to be given 

by 




1/3 




(6.7) 


where qs is a constant. In this case of discontinuous solutions, the jump conditions relate 
the dependent variables at q = q^ to those at q = q~ and are given by 
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( 6 . 8 ) 
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Key locations in the similarity solutions are the points at which the dependent variables 
transition from the new steady state to an unsteady pulse and then from this unsteady 
pulse to the original steady state. The family of slowest moving characteristics associated 
with the new source is given by 

i = M(5-ys(sTT)) = ™(s_ys(sVTj). 

The transition between the steady and unsteady portions of the solution must occur at 
a singular point of D to permit the gradient of the solution to change discontinuously. 
Thus in this case q/fh = A{S— S {S — 1) j /3 and so we find that 


dz 3z 
dt 4< ’ 


( 6 . 10 ) 


This implies that the transition point occurs at a constant value of the similarity variables, 
V = Vci, given by 

ryei = 1{S- VSiS-1)) . (6.11) 

Likewise the family of fastest moving characteristics is given by 

di = ^(s+yS(SVT)) = ™(s + VS(SVl)). (6.12) 


and so the boundary between the unsteady pulse and the steady solution associated with 
the original buoyancy flux occurs at a constant value of the similarity variable, rj = r]c 2 
given by 

Vc2 = {s + VS{S-l)) , (6.13) 

provided the motion is due to a decrease in the source strength {P > 1). If the source 
strength increases then the characteristics associated with the new release overtake those 
due to the original source and, as we show below, the motion forms ‘shocks’. 


6.1. Increase in buoyancy flux 

These flows correspond to a new scaled buoyancy flux such that < 1. For rj < rjd, the 
similarity solution is given by 



_ fl 
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fj2- 

’ v^J 


(6.14) 


where r]ci is given by (16.111) and corresponds to the slowest moving characteristics associ¬ 
ated with the new buoyancy flux. The leading edge of the unsteady solution corresponds 
to a shock at rj = rjs such that, for rj > rjs, the solution is given by the steady state 
associated with the original buoyancy flux given by (12.121) . 

Constructing the similarity solution then requires the computation of the solution for 
Vci < V < Vs, where Vs is to be determined as part of that solution. We note that this 
domain includes a location where q/fh = 4/3. Here there is the potential for a con¬ 
tact discontinuity where the reduced gravity is discontinuous or even unbounded, but 
the volume and momentum fluxes are continuous. The numerical problem is therefore 
a boundary-value, ordinary differential equation with an internal critical point, which 
we solve using a numerical shooting technique. Our method of solution proceeds as fol¬ 
lows. First numerically integrate from r] = ijci using the local series expansion given in 
appendix o This expansion provides a solution local to the critical point and entails 
an undetermined constant, Ca- Given a value of Ca, the numerical integration can be 
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Figure 6. The scaled width, B{z) = the scaled velocity, W{z) = z 

and the buoyancy flux, F as functions of the distance from source, z at dimensionless times 
t = 2.4,4.9, 7.4 and 9.9 for an increase in buoyancy flux {F = 0.05) at t = 0 and with shape 
factor S — 1.1. In (b) the scaled velocity, W = , corresponding to the steady velocity field 

of the original buoyancy flux at the source, and the scaled velocity W = , correspond¬ 

ing to the steady velocity field associated with the new buoyancy flux, are plotted in dashed 
lines. In (c) the original buoyancy flux, F = Fq — 0.05, and the new steady state buoyancy flux, 
F — Fi = 1, are plotted in dashed lines. 


continued until rj = 77^1 at which q/fh = 4/3. The solution is also numerically integrated 
from the shock at the leading edge. Given a shock location ps, the jump conditions (16. 8|) 
provide the conditions at r] = rjJ and then the solutions may be numerically integrated 
until r] = r]rn 2 at which q/rh = 4/3. The constant Ca and the shock location r]s are then 
iteratively adjusted until ijmi = r]m 2 and Q{r]mi) = Q(jim 2 ) (noting that this ensures that 
the momentum flux is continuous at this location as well), and when these conditions are 
satisfied, this provides the entire solution. 

The numerical integration of the governing partial differential equations is plotted in 
figure |6] at various instances of time and the underlying similarity form of solution is 
plotted in figure [7l We note that there is excellent correspondence between the two, 
as evidenced by the overlap between the similarity solution and the direct numerical 
computations (the curves in figure Hare virtually indistinguishable). We observe in figure 
M that an increase in buoyancy flux at the source leads to a broadening of the width of 
the plume in the unsteady pulse before the new steady state is established. Notably the 
velocity field is increased relative to the flow associated with the original buoyancy flux. 
Together these lead to the surprising variation in the buoyancy flux as it changes from the 

















24 M.J. Woodhouse, J.C. Phillips & A. J. Hogg 



Figure 7. The similarity solution for the volume flux, q{g), the momentum flux, rh{g), and the 
buoyancy flux, f{g) as functions of the similarity variable g for shape factor S' = 1.1 and an 
increase in the source buoyancy flux T = 0.05 (plotted in solid lines). Also plotted are results 
from the direct numerical integration of the governing equations (dashed lines), although the 
two sets of curves are so close in values that they are virtually indistinguishable. The values 
gd, gm and gs are also marked. These corresponds, respectively, to the boundary between the 
time-dependent part of the solution and the steady state associated with the new buoyancy flux 
at the source, the location of a contact discontinuity and the location of a shock that marks 
the interface between the unsteady evolution and the steady steady associated with the original 
buoyancy flux at the source. 

original value oi F = 0.05 to the new value of = 1; the variation is not monotonic and 
within the unsteady pulse the buoyancy flux overshoots its new value before subsequently 
decreasing to attain it (see figureE]:). The similarity structure is plotted in figure^ where 
the similarity variables are plotted between the leading and trailing edges of the unsteady 
region. These correspond to a shock at g = ps and a point where the gradient changes 
discontinuously at g = gd- In between there is a contact discontinuity at 77 = 77 ^ where 
the volume and momentum fluxes are continuous, and the buoyancy flux is unbounded. 

6.2. Decrease in buoyancy flux 

When the buoyancy flux decreases relatively weakly {F < Fm = 6.9), we construct the 
similarity solution between the two critical points, gd and gc2, given by (16.111) and (I6.13|) . 
respectively. For g > gc2, the solution corresponds to the steady state associated with the 
original buoyancy flux, whereas for g < gd, the solution corresponds to the steady state 
associated with the new buoyancy flux. As described above, the boundaries between the 
steady states and this unsteady pulse are characteristics that propagate at the fastest 
and slowest rates. At some point, g^, (jid < g-m < gc2) the similarity solution reaches a 
state in which q/fh = 4/3 and there is a contact discontinuity. 

We construct the solutions as follows: we integrate from g = gd, initiating the numer¬ 
ical solver with the local expansion derived in appendix [C] which entrails an adjustable 
constant, Cai- We numerically integrate until g = gmi at which point q/m = 4/3. We 
also numerical integrate from g = gc2, initiating the solution using a local expansion de¬ 
rived in appendix [cl which features another adjustable constant Ca 2 - We integrate until 
g = gm2 at which point q/fh = 4/3. The solution is then found by iteratively adjusting 
Cal and Ca2 until gmi = gra2 and q{,gmi) = q{gm2)- 
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Figure 8. The scaled width, B{z) = zQ/'/m, the scaled velocity, W{z) = z~^l^MIQ 
and the buoyancy flux, F as functions of the distance from source, 2 , at dimensionless times 
t = 4.5, 9.5, 14.5 and 19.5 for a decrease in buoyancy flux [T = 2) aX t = Q and with shape factor 
S = 1.1. In (b) the scaled velocity, W = 2“^^®, corresponding to the steady velocity field of the 
original buoyancy flux at the source, and the scaled velocity W = , corresponding to 

the steady velocity field associated with the new buoyancy flux, are plotted in dashed lines. 


The numerical results from direct integration of the governing partial differential equa¬ 
tions are plotted at various instances of time in figure |8l Here we observe that the volume 
and momentum fluxes evolve continuously, in contrast to the evolution following an in¬ 
crease in source strength ( 116.II) . During the unsteady evolution from the original state to 
the new one, the radius of the plume decreases and the velocity increases. The buoyancy 
flux, however, does not monotonically vary from the new values {F = 1) to its original 
value {F = 2). Instead it initially decreases (see figure |5J:). The similarity solution for 
the unsteady pulse is plotted in figure [9] between the leading and trailing characteristics 
(FcI < r] < r]c2)- This solution features a contact discontinuity at ij = ijm, although the 
mass and momentum fluxes remain continuous at this point. In figure IHl we have plotted 
both the similarity solution and the rescaled results from the direct numerical integra¬ 
tion of the governing equations, and we note that the curves for each of the fields are 
indistinguishable in this plot. 

The morphology of the solution changes for larger decreases in the buoyancy flux. 
When F > Fm, we find that the similarity solution features another critical point at 
r] = r]c3 at which det(D) vanishes and the system potentially becomes singular. Local 
analysis of the behaviour close to this new critical point indicates that non-integer powers 
in the series expansion are not possible here; the determined power a is negative and so 
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Figure 9. The similarity solution for the volume flux, q{g), the momentum flux fh(ri) and the 
buoyancy flux, /(ry) as functions of the similarity variable g for shape factor S = 1.1 and a 
decrease in the source buoyancy flux T = 2 (plotted in solid lines). Also plotted are results 
from the direct numerical integration of the governing equations (dashed lines), although the 
two sets of curves are so close in values that they are virtually indistinguishable. The values 
gd, g-m and ? 7 c 2 are also marked. These correspond, respectively, to the boundary between the 
time-dependent part of the solution and the steady state associated with the new buoyancy flux 
at the source, the location of a contact discontinuity and the location of a the interface between 
the unsteady evolution and the steady steady associated with the original buoyancy flux at the 
source. 

in order to ensure the fields are bounded, we must enforce Caz = 0. This implies that the 
dependent variables pass smoothly through this critical point. However, having attained 
a state in which q/fh > 4(5 -|- \JS [S — l))/3, the only way to connect to the rest of the 
solution is via an internal shock at g = ps- Our method for constructing the solution 
then proceeds as follows. We integrate numerically from the critical point at 77 = gci, 
initiating the solution using the local series expansion about this point and introducing an 
adjustment constant, Cai, to a location g = gmi at which there is a contact discontinuity 
[q/fh = 4/3). We also numerically integrate from g = gc 2 , initiating the solution using 
the series expansion with constant Ca 2 - This constant is adjusted so that an internal 
critical point is reached at g = gcs {gcs < gc 2 ), where we enforce the solvability condition 
given in ini. We may then integrate further; we smoothly pass through the critical 
point and insert a shock at 77 = rys < gc 3 and then integrate until q/fh = A/'i at g = 77 ^ 2 - 
This leaves two adjustable constants, namely Cai and gs^ which are iteratively adjusted 
until gmi = gm 2 and q{gmi) = q{gm 2 ) to give the complete solution. 

The results from the numerical integration of the governing equations are plotted in 
figure [To] for a large decrease in buoyancy flux at various instances of time. We note that, 
as with the weaker decreases in flux, the plume responds by narrowing and accelerating 
in order to adjust back to its original state. However a small internal shock is also 
developed. For these parameter values {P = 20 and S = 1.1), the shock is of relatively 
small magnitude and it generates discontinuities in the velocity and width fields, with 
the reduced gravity remaining continuous. From this figure it is not possible to observe 
the presence of the internal critical point (77 = gcs) because all of the variables and their 
derivatives are continuous. However, it can be confirmed that there is an internal region 
within which q/fh > 4(5 -I- S {S — l))/3. The associated similarity solution within the 
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Figure 10. The scaled width, B{z) = zQ/\/m, the scaled velocity, W{z) = z~^^^ MIQ 
and the buoyancy flux, F as function of the distance from source, z, at dimensionless times 
t = 4.5, 9.5,14.5 and 19.5 for a decrease in buoyancy flux F — 0.05 at t = 0 and with shape 
factor S = 1.1. In (b) the scaled velocity, W = corresponding to the steady velocity 

field of the original buoyancy flux at the source, and the scaled velocity W = , 

corresponding to the steady velocity field associated with the new buoyancy flux, are plotted in 
dashed lines. 

unsteady pulse is plotted in figure [TTJ As with the weaker decrease in buoyancy flux, 
this plot features continuous transitions at the leading and trailing edges (77 = ijd and 
V = ^c 2 , respectively) and a contact discontinuity [rj = rjm). Additionally, there is an 
internal critical point (j] = rjcs) and an internal shock (77 = rjs). The rescaled results 
from the numerical integration of the governing equations are overlain on the similarity 
solutions and again they are virtually indistinguishable, confirming the presence of this 
similarity solution in the unsteady dynamics of the plume. 


7. Discussion and conclusion 

Steady plume models are applicable on time scales that are long compared to the eddy 
turnover time that characterises transient turbule nt features, and on time scales shorter 
than the characteristic time for source variations ( Scase et a/.l 120066 b In many applica¬ 
tions, the effect of variations of the source fluxes of mass, momentum or buoyancy on 
the plume dynamics are of fundamental importance. Here, unsteady models are essential. 
For an example, during volcanic eruptions the source conditions can fluctuate in time 
due to unsteadiness in the physical processes occurring in the volcanic conduit and at the 
vent. End members of the source unsteadiness are the initiation phase (when an eruption 
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Figure 11. The similarity solution for the volume flux, q(rf), the momentum flux rh{rf) and 
the buoyancy flux, f{g) as functions of the similarity variable g for shape factor S' = 1.1 and a 
decrease in the source buoyancy flux = 20 (plotted in solid lines). Also plotted are results from 
the direct numerical integration of the governing equations (dashed lines), although the two sets 
of curves are so close in values that they are virtually indistinguishable. The values g^i, g-m, 
gs, gc 3 and gc 2 are also marked. These corresponds, respectively, to the boundary between the 
time-dependent part of the solution and the steady state associated with the new buoyancy flux 
at the source, the location of a contact discontinuity, the location of an internal shock between 
flow states, the location of a continuous transition and the location of a the interface between 
the unsteady evolution and the steady steady associated with the original buoyancy flux at the 
source. 


first produces material that becomes buoyant and ascents through the atmosphere) and 
waning phase (when the eruption comes to an end either through a gradual reduction 
in the rate at which material is erupted or an abrupt cessation of the activity) of the 
eruption. 

The unsteady model of turbulent plumes analysed by Scase et ^ (20063) and IScas3 
(20 o3 ) identified a key feature of the unsteady response of plumes to abrupt changes 
in the source buoyancy flux; the plume adjusts to the new source c onditions through a 
‘pulse’ that propagates through the plume. However, the model of Scase et ^ ( 200661) 
adopts top -hat profiles for t he me an axial velocity and was shown t o be ill-posed in the 
ana lysis o f Scase fc HewittI ( 2012h . with the numerical solutions of Scase et al. ( 20066 1 
and Scase ( 2009ll attained only due to significant numerical visco sity. _ 

The diffusive regularisation of the unsteady model proposed bv IScase fc Hewitt (20121 
has an intuitive physical basis: turbulent eddies have a vertical l ength scale over which 
different levels of the plume are connected. Numerical solutions (jScase fc Hewittll2012ll 
suggest t hat the phenonr e nolog ical model of the turbulent diffusion of momentum advo¬ 
cated by Scase fc Hewitd ( 201211 results in a well-posed unsteady model, but the analysis 
presented here shows that this modification of the governing system of equations ren¬ 
ders the classical power-law solution for steady plumes unstable. As the power-law solu¬ 
tions have strong empir ical support (e.g. Morton et allllOS^: IPapanicolaou fc List 1988: 
Shabbir fc GeorrallOOdll . we conclude that the diffusive regularisation of IScase fc HewittI 


(I2OI2I) is not appropriate. 

Our analysis identifies a different physical process, a difference in the rates of transport 
of the cross-sectionally averaged mass and momentum of the plume (referred to as type 
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I dispersion by Craske fc van ReeuwiikI 2015311) that occurs due to non-uniform radial 
profiles of the mean axial velocity in the plume, and we have shown that including this 
process through a momentum shape factor leads to a well-posed system of equations. 
However, this requires the top-hat description of the radial profiles for the axial velocity 
in the plume, which has been applied extensively (for convenience) in models of steady 
plumes, to be replaced with a description of the radial variation. The resulting integral 
model for unsteady plumes introduces only one add itional parameter, the momentum 
shape factor, over the classical steady plume model (iMorton et al.lll9^ . Furthermore, 
the (ensemble averaged) radius of the plume remains dependent on the entrainment 
coefficient alone, so the unsteady model does not preclude calibration of the entrain¬ 
ment coefficient from laboratory experiments that measure the steady plume width (or 
spreading angle). To determine the mome ntum shape factor, the radial profile of the ax¬ 
ial ve locity can be measured directly (e.g. Papanicolaou fc Li3]ll988l : IShabbir fc Georg3 
I 1994 II and the shape factor computed from (1^ . 

Explicitly including a momentum shape factor that differs from unity results in a 
strictly hyperbolic system of eq uations that gover ns the plume dynamics. Th is is appeal¬ 
ing a s laboratory experiments ( Scase et al. 2008^ and numerical modelling ( Scase et al 


I 2 OO 9 II show that the adjustment of the plume to changes in the source conditions occurs 
through the propagation of an unsteady pulse, as predicted by our non-diffusive hyper¬ 
bolic model. The development of the unsteady pulse following an abrupt change in the 
source buoyancy flux is described by a similarity solution of the hyperbolic system of 
equations. The construction of the similarity solution allows us to identify three quali¬ 
tatively different regimes of the unsteady evolution. Following an increase in the source 
buoyancy flux, the pulse takes the form of a localised increase in the plume width with 
a leading discontinuity. If the source buoyancy flux is reduced then the plume width 
narrows. For a relatively strong reduction in the source buoyancy flux, an internal shock 
occurs in the similarity solution, whereas no internal shock is found in the source buoy¬ 
ancy flux is reduced by a smaller amount. We expect diffusive processes will act to smooth 
locally the sharp gradients that appear in solutions of the hyperbolic model. However, 
hyperbol ic models have been shown to capture the dominant flow dynamics in many 
se ttings dWhitham 1974 ). 

ICraske fc van Reeuwiik (20153) identify, from direct numerical simulations, the disper¬ 
sion of momentum as a fundamental feature of turbulent jets and construct an integral 
model to describe unsteady j ets that includes a description of the non-uniform radial 
profile of the vertical velocity ( Craske fc v an Reeuwiikll20153) through shape factors. In 
the model of ICraske fc van ReeuwiikI (|20156ll the integral conservation equations are de¬ 
rived from t he point-wise momentu m and energy conservation equations, following the 
approach of Priestley fc Ball! (1955), and an energy shape factor rather than a momen¬ 
tum shape factor is introduced. The resulting system of integral conservation equations 
shares some features with the unsteady integral model proposed here (albeit for a jet 
rather than a plume), in particular the hyperbolic structure of the system of equations. 

For jets, the momentum of the flow as it is ejected from a source drives the mo¬ 
tion, and there is significant ev olution of the radial profile of the axial velocity (re¬ 
ferred to as type H dispersion bv ICraske fc van' Reeuwii3l2015 30) as the flow develops 
away from the source. The evolution of the shape of the radial profiles of axial veloc- 
ity and buoyancy have bee n incorporated i nto a non-constant entrainment coefficient by 


Kaminsk i et al. (l2005h and Carazzo_el_a|. ( 20061) . Furthermore, in the unsteady integral 
model of ICraske &: van ReeuwiikI (|20156l) . the deviation of the radial profiles from self¬ 
similar forms gives rise to diffusive terms in the system of integral conservation equations. 
Our analysis shows that a description of momentum dispersion, through a shape factor 
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that differs from unity, is sufficient to obtain a well-posed model of unsteady plumes. 
Therefore, while for jets it is necessary to account for the evolution of the radial profile 
of axial v elocity to self-s i milar form, for unsteady pure plumes the classical entrainment 
closure of Morton et al. ( IQSd) remains applicable. 

The mathematical model we present allows predictions of unsteady plume dynamics 
to be made and, while our study has focussed on changes to the source buoyancy flux, 
the effect of general temporal variations in the source conditions can be examined. Fur¬ 
thermore, our framework can be applied to describe the unsteady dynamics of plumes 
in industrial and environmental settings, where additional physical processes such as an 
ambient flow, thermodynamics and particle transport have a strong influence on the evo¬ 
lution of the plume. 


We are grateful to Professor R.R. Kerswell and Dr C.G. Johnson for valuable discus¬ 
sions. This project has received funding from the European Union’s Seventh Programme 
for research, technological development and demonstration under grant agreement No. 
308377 and from the NERC through grant NE/I01554X/1. 


Appendix A. Conservation of energy for unsteady plumes 

In the derivation of our integral model for unsteady plumes, we made direct use of the 
point-wise continuity equation t o obtain an in t egral equation for conservation of mass, 
following the approach used bv iMorton et al. ( 1956h for steady plumes. This requires 
a priori a choice to be made for a representative plume radius, &, and an entrainment 
closure to describe the mixing of ambient fluid into the plume. 

An alterna tive a pproa c h, pioneered by Priestley fc dlQS Sil for stead y flows and 
developed bv IFoxI (|i 970I1 . Kaminski et al. ( 2005 ) and Carazzo et al\ ( 20061) . is to sub¬ 
stitute conservation of axial kinetic energy for the continuity equation in the set of 
point-wise conservation equations. As the equation for conservation of axial kinetic en¬ 
ergy is derived as a (axial velocity) moment of the equation for the conservation of axial 
momentum, rewritten in conservation form through application of the continuity equa¬ 
tion, taking the conservation of axial kinetic energy together with the conservation of 
axial momentum provides no additional information over the set of point- wise conserva¬ 


tion equations with continuity and conservation of axial momentum (see lMortonlll971 


for a detailed discussion). The integral expression for conservation of axial momentum 
and conservation of axial kinetic energy can be manipulated into a form that expresses 
conservation of mass in the integral sense and, from this, an expression for turbulent 
entra i nment is obtained (IPriestlev fc Balllll955 : FoxI 1970l: Morton 1971 : Kaminski et al 


20051: Carazzo et al. 2006 ). However, the integral equations that result from the use of 
the conservation of axial kinetic energy differ from those obtained when the continuity 
equation is used directly. Indeed, the differences occur due to the application of turbu¬ 
lent closures in different terms following t he cross-sectional int egration of the point-wise 
conservation equations: in the approach of iMorton et al. 1 19561) the turbulence closure is 
applied as a inflow velocity in the kinematic condition applied at the plume boundary, 
whereas the closure is applied to the turbulent fluctuation te rms that occur in the co n- 
servation of axial moment um equa t ion w hen the approach of Priestley fc Balll ( 1955[) is 
used. Finally, as noted bv iMortonI (|l97l[) . both closures introduce parameters that can 
only be det ermined empirically. 

Recently, ICraske fc van Reeuwiik ( 20150116 ) used the conservation of axial kinetic en¬ 
ergy is place of the continuity equations in the set of point-wise conservation equations in 
their model of unsteady jets. In this appendix we adopt a similar approach for buoyant 
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plumes. We note that our approach differs as we again specify the plume radius r = b{z, t) 
as a measurable surface in the plume, although we do not define this precisely, and in¬ 
troduce shape factors to account for non-uniform radial profiles of mean flow quantities. 
Thus, we define the volume flux, axial momentum flux, and the flux of mean axial kinetic 
energy as 


b^W = 2 rwAr, Sh^W'^ = 2 


ruP' dr, and SePW^ = 2 




dr. 


(Ala, b, c) 

respectively, where, in (lA Ib l. Se is an energy flux shape factor. 

Integration of the point-wise Reynolds-a veraged conservation equations for axial mo¬ 
mentum (12.161) , axial kinetic energy (see ICraske fc van ReeuwiikI l2015all and reduced 
gravity (I2.1cl) gives the following set of integral equations. 


4 ^ (Sb^W^) = PC + Mp ), (A 2a) 

(Sb^W^) + ^ (Seb^W^) = 2(t>b^CW 
ot ^ ^ oz ^ ^ 

g 

+ {Pm + Pf + Pp) ~ {Ef + Ep) , (A 26) 

+ ,A2c) 

where, following ICraske fc van ReeuwiikI (|2015ai r^. the correlation fluctuation terms and 
non-hydrostatic pressure terms, that are formally higher order in the plume slenderness, 
R/E[ , than the mean fluxes, are retained and given by 


Mf =2 rw''^ dr, Mp = 2 rpdr, Ef = A dr, 

7n Jn Jn OZ 


Ep=A 


Pp =4 


rpw dr, Pm = 4 / ru'w‘ 


■ dw 
dr 


rdw 


P _dw ^ 

rp—dr, Bf=2 


P _ 

/ rg'^w' dr. 
JQ 


(A3) 


Here p denotes the averaged deviation of the pressure in the plume from hydrostatic 
pressure. Note, in (lA 21) we have neglected terms in which flow quantities are evaluated 
on the plume boundary. 

Manipulation of the integral equations for conservation of axial momentum (I A 2al) and 
conservation of axial kinetic energy (lA 26^ leads to an expression akin to conservation of 
mass, 




(A 4) 


where the entrainment rate ke is given by 


ke = 



bC 

1 1 1 


5^ 

\ ^ 

V ) 


6 W2 ' 

U 

SeJ 

' dz 


{sPw^ 


05' 52 0 /S, 


2SeW 06 2Se dz V 52 


Pm + Pm+Pf S ^ f l\4- , ]\/r \ I ^ ^ I'JZ I 


(A 5) 


The form of the mass conservation equation (1X41) is similar to the expression pro- 
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posed by Craske fc van R.eeiiwiik ( 2015a but differs due to different choices for the 
plume w idth. The entrainment rate has a similar form to that of Craske fc van R eeuwiikI 
( 2015alf^ for jets, although in (lA 51) the first term on the right-hand-side does not appear 
for jets where G" = 0. A similar term, which results in an entrainment rate that depends 
on the local Richardson number of the plume (given by Ri = hG o c curs i n the 

entrainment rate in the models that adopt the approach of Priestley_fc Bah ( 1955[l (e.g. 
Priestley fc Ball 1955 : FoxIIlQTOl: Kaminski et al.ll2005t ICarazzo et aU 2006ll and has been 
shown t o be important when d e scribing the near source development of bu oyant jets into 
plumes ( Kaminski et al. 2005 : Carazzo et al. 2006t Ezzamel et al.l[2015ll . Further, the 
entrainment rate here includes a contribution from temporal ch anges in the momentum 
shape f actor that does not appear in the entrainment rate of Craske fc van ReeuwiikI 

(l2015airai . 

For turbulent plumes from pure plume sources at distances sufficiently far from the 
source, experiments suggest the radial profiles of mean axial velocity and reduced grav- 
ity, and the second-o rder turbulent velocity statistics, have reached a self-similar form 
( Ezzamel et all 2015ll. The shape factor s can th en be taken as constant values. Further¬ 
more, following Craske fc van ReeuwiikI ( 201561) and for simplicity, we take the turbulent 
production and transport terms and non-hydrostatic pressure terms to be such that the 
entrainment rate is given by. 


ko -(- 


S-: 


bC 


bW^ 


1 . 


dz 


(Sb^W^) 


(A 6) 


with fco constant. 

Steady solutions of the plume model with the non-constant entrainment rate 
with pure plume boundary conditions (Q(0) = 0, M(0) = 0, F{0) = Fq) are given by. 


Q = 


6ko f9koV^^ fFoV^^ \ 


5 

M = 

F = Fo. 


9fco 

10 


10 y V -S' 

2/3 /p \ 2/3 r 


1 - 


8S 




1-4/3 


.5/3 


1 - 


85' 

55« 


(5 - 52 - y. + Se) 


1 -2/3 


. 4/3 


{A 7 a) 

(A 76) 
(A 7c) 


Therefore, the non-constant entrainment coefficient does not change the spatial variation 
of t he steady solu t ions, although the pre-factors in the power-law forms are changed from 
the Morton et al. ( 195(tI ) solutions, unless the shape factors take values such that 5^ — 5— 
Se + cj) = 9- We note that the factors could diverge if 1 — 85 (5 — 5^ — -I- 5e) /55e = 0. 

The radius of the plume for the non-constant entrainment coefficient is given by 


6 = 


6fco 

5 




(AS) 


and so the plume radius increases linearly with height as in the model of Morton et al\ 
i 19561) . However, in contrast to the model with a constant entrainment coefficient, the 
spreading angle of the plume depends on the the shape factors for momentum, buoyancy 
and energy (unless — S — Se + (j) = 0)- 

An alternative approach is to augment the system of integral equations (lA 2p with the 
integral expression for conservation of mass (I2.10ap . We can then interpret the integral 
expression for conservation of axial kinetic energy as governing the evolution of the mo¬ 
mentum shape factor and, in order to fully specify this evolution, a closure is required 
to describe the evolution of the energy shape factor. A higher velocity moment of the 
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point-wise conservation equations for mass and momentum would provide an integral 
equation governing the evolution of the energy shape factor, but would, of course, intro¬ 
duce a new shape factor. Therefore, a turbulence closure must be i nvoked at some level 
within the hierarchy of equations. Laborator y experiments (e.g. iPapanicolaou fc List 
1988t Wang fc Law! 2002 : Ezzamel et aLllioTBh suggest the momentum shape factor can 


be taken as a constant value for fully developed turbulent plumes sufficiently far from 
the source, and this represents a turbulence closure at the level of conservation of mass, 
momentum and buoyancy. 


Appendix B. A phase plane analysis of steady solutions of the plume 
model with diffusive terms. 

In ^ the steady solutions of the integral model for unsteady plumes with phenomeno¬ 
logical diffusive terms included to describe turbulent mixing by eddy diffusion were shown 
to be linearly unstable to small perturbations. Therefore, while steady solutions of the 
diffusive system of equations can be found analytically ( s ee 13.31) and are modifications of 
the well-known steady solutions given bv iMorton et al\ (1195611 . the solutions cannot be 


realised. However, the linear stability analysis in IjS] does not investigate other possible 
steady states. In particular, it could be possible that other steady solutions that do not 
enforce pure plume boundary conditions exist, and that these could be attracting states 
for the steady diffusive system of equations. 

Here we examine the steady plume equations with diffusive terms by treating the 
system of equations as a dynamical system. To this end, we allow the volume flux Q to 
be the independent variable in the system rather than the distance from the source. We 
note that, in an unstratified ambient the buoyancy flux is conserved, and so F = Fq with 
Fo a specified constant. 

We consider first the diffusive term proposed by IScase fc Hewitt (2012 h (as given in 
the unsteady momentum conservation equation 13.2|) . The equation for the conservation 
of mass in the steady state allows us to write 


and therefore, from the conservation of momentum we find 

+ ,91^" + (3, + S)0^ + 2.M --99^ 

^ dQ2 dQ 2fcM3/2’ 


(HI) 


(B2) 


here treating M as a function of Q. Equation 
introducing new dependent variables 


can be further manipulated by 


X = Q-^/^M{Q) and Y = ^ M{Q), 

dQ 5 


(B3) 


and an independent variable ^ with Q = exp (^) which allows the system to be written 
as the following autonomous system, 


dX 


= T, 


^ V + f ^ V =-^x-3/2. 

d^ X V 5 J \25 5 J 2k 


(B4a) 

(B46) 


The coupled nonlinear system 


has a single fixed point at T = 0 and X = Xq 
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Figure 12. Trajectories i n the phase-plane of sol utions to the steady plume equations with the 
diffusive regularization of IScase fc HewitH l|2012l l with fb = 1, A: = 0.1, k = 0.1 and S = 1. 
The autonomous system has a single fixed point ( denoted by the black point on the figure) that 
corresponds to the steady pure plume solution of iMorton et al\ (Il95fil ~) with a modification for 
the diffusive term (given bv l3.3ll . The fixed point is a saddle, with two stable directions and two 
unstable directions, shown as bold solid lines on the figure with arrows denoting the stability of 
the trajectories. Additional trajectories, corresponding to solutions with non-pure plume initial 
conditions, are shown as thin grey lines. 


with 

which, after returning to the original variables, is the steady solution given by (13.31) . 
However, other solutions of the autonomous system for arbitrary initial conditions can 
be readily found by numerically integrating the system (El. The trajectories in the 
phase plane corresponding to solutions of (lB4l) with a momentum shape factor S' = 1 
are shown in figure[T2l We see that the fixed point is a saddle with two stable trajectories 
(such that initial conditions precisely on these trajectories converge to the fixed point) 
and two unstable trajectories. For initial conditions that are perturbed away from the 
pure-plume conditions (i.e. Q{z = 0) = M(z = 0) = 0) at the fixed point, the trajectories 
move the solution away from the fixed point and we find either M{Q) —?> 0 or M (Q) —>■ oo 
along the phase-plane trajectories. Taking a momentum shape factor that differs from 
unity does not change the topology of the trajectories in the phase plane (figure [H. 

A local analysis of trajectories that are perturbed away from the fixed point reproduces 
the results of the linear stability analysis dSH). Indeed, taking X = Xq + Ai(^) and 
Y = yi(f) with I All <C 1 and |Yi| <C 1, we find 

Ai\ / 0 1 \ /Ai 

YiJ \1/5 + 2 S/k A/S, + S/kJ\Yi 
and the eigenvalues of the Jacobian matrix are 

5S + 4k± V25S'2 -b 240 S'k - 




o± = 


10k 


(B7) 
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Figure 13. The stable and unstable trajectories through fixed points in th e phase-plane of 
solutio ns to the steady plume equations with the diffusive regularization of IScase fc HewittI 
(I2OI2I ') with Fq = 1, k = 0.1 and k = 0.1 for momentum shape factor S = I (solid lines), S = 1.5 
(dashed line) and S = 2.0 (dotted line). 


The corresponding eigenvectors are 


v± = (l,a±) 


(B8) 


Therefore we find Xi ^ expa^. 

The analysis presented above can be applied to the alternative form of the eddy dif¬ 
fusion term in the momentum balance given in (13.81) for unsteady conditions. We again 
find a single fixed point that corresponds to the steady solution (1^ . The topology of 
trajec tories in the phase plane are qualitatively similar to those of the IScase fc Hewitd 
( 2 OI 2 I) diffusive term (figures [TTl and [T51) with the single fixed point being a saddle and 
trajectories from arbitrary initial conditions have either M(Q) —>■ 0 or M{Q) —>■ 00 as 
Q —>■ 00 . Taking a momentum shape factor that differs from unity does not change the 
topology of the phase portrait (figure fT^. 

The conclusion of this analysis is that, for each of the proposed axial diffusion terms 
()3.1l and l3.7l) . the steady solutions for pure plume boundary conditions are not stable and, 
furthermore, for more general boundary conditions, the steady solutions do not evolve 
towards the pure plume solution in the far-field. This is in contrast to steady solutions 
of the non-diffusive system wit h general boundary conditions that con verge to the pure 
plume solution in the far-field ( Morton 19591 : Hunt fc Kave 2001 . 2005h . 


Appendix C. Local series expansions at critical points 

In this appendix we examine similarity solutions governed by 

Dv^ = b, (C 1) 

ar] 

where D and b are the matrix and vector defined in (16.31) . We examine the solutions local 
to a critical point, rj = rjd (i = 1, 2), of this differential system, where det(D) vanishes. 
It is convenient to write v = \og{r]/r]ci) and to write the following expansion of the 
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Figure 14. Trajectories in the phase-plane of solutions to the steady plume equations with the 
diffusive regularization (|3.8|) with fb = 1, fc = 0.1, k = 0.1 and 5 = 1. The autonomous system 
has a single fixed point (denoted by the blac k point on the figure) that corresponds to the steady 
pure plume solution of iMorton et al\ lIlOhfiTl with a modification for the diffusive term (given by 
The fixed point is a saddle, with two stable directions and two unstable directions, shown 
as bold solid lines on the figure with arrows denoting the stability of the trajectories. Additional 
trajectories, corresponding to solutions with non-pure plume initial conditions, are shown as 
thin grey lines. 



Figure 15. The stable and unstable trajectories through fixed points in the phase-plane of 
solutions to the steady plume equations with the diffusive regularization (1^ with Fb = 1, 
k = 0.1 and k = 0.1 for momentum shape factor S' = 1 (solid lines), S = 1.5 (dashed line) and 
S = 2.0 (dotted line). 

dependent variables, the matrix D and the forcing vector, b, 

*2"^ rv—I—1 

9 — Qo + + U q2 -b U Qa -b U Qa-l-l + • • • 5 

D = Dq -b u Di -b Dq, -b ..., 
b — bo P vb\ -b v°‘bci -b .. •, 


(C2a) 

(C 26 ) 

(C2c) 
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where a is not an integer. This form of expansion series is possible because det(Do) = 0. 
We substitute the series (1C 21) into m and in the regime <C 1, balance terms in 
powers of ly. At 0(1) and at 0(z^““^) we find that 


DoQi = bo and aDoq„ = 0. (C3a, 6) 

The matrix Dq is singular; thus we can find vectors e and such that DqB = 0 and 
e^Do = 0, respectively. This means that from (1C 3b ') we deduce the solvability criterion 

e^bo = 0, (C4) 

and that = Cie + q^ and q^ = Cae, where q^ is a particular solution of (1C 3h l and 
Cl and Ca are constants. From a balance of terms at 0{iy) and 0{v°‘) we deduce 

Diqi — bi and {Datli + Diaq^) = e^ba. (C 5a, b) 

The first of these equations (1C 5b l determines the value of the constant Ci, while the 
second (jC 5b l determines the non-integer power a. 

We may apply this formulation to any of the critical points, but the locations that 
are of most interest are the boundaries between the steady and unsteady portions of the 
solution. First we analyse the local behaviour near to qc 2 = (^S + S {S — 1)^ /3, 

which corresponds to the boundary between the original source and the unsteady pulse 
within the similarity solution. At this point the solution is given by 

90 = (/3o,/3o./3o) . where /3o = ^ (*5 + x/S" (5 - 1)) . (C6) 

The vectors e and are given by 

= (^fS (2S-1 + 2^S{S-1)) , I (25 - 1 + 2v'5(5-l)) , l) , (C 7a) 
= (-1 (25 - 1 - 2v'S {S - 1)) , 1 , 0 ). (C76) 

The particular solution is given by 

g/ = (-/3o,-2/3o^-3/3o'), (C8) 

the constant Ci = 0 and the exponent, a is 

(C9) 

8 16v'S'(S'- 1) 

The local expansion is bounded as —>■ 0 provided a > 0, a condition that demands 
S > 25/24. We note that this condition is identical to the criterion for linear stability. Its 
origin is identical; it comes from the requirement that the unsteady adjustment remains 
bounded at its leading edge. 

The solution at the trailing edge of the unsteady portion of the solution is broadly 
similar. Expanding the dependent variables close to rjd = 4 S {S — 1)^ /3, we 

find that the solution is given by 

90 = where /3i = | (^- (5 - 1))(CIO) 

The particular solution is given by 

9/ = (~/3i,—3/3o) , 


(Cll) 
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the constant Ci = 0 and the exponent, a is 

13 5(25'-!) 

a =-^- \ ’ ■ 

8 ley's” (S'- 1) 


(C12) 


Appendix D. Separable similarity solution 

IScase et~^ (12006 6li identified a ‘separable’ similarity solution that emerged as an exact 
solution to their time-dependent governing equations and, although it did not satisfy the 
boundary conditions exactly, it appeared to play a role in ‘organising’ the underlying 
dynamics when the buoyancy flux at the origin was abruptly reduced by a relatively 
large factor. In this appendix we derive the analogue of their solution in our regularised 
system of governing equations and discuss its relevance for the solutions that arise when 
the source buoyancy flux is abruptly changed. 

The separable solution emerges as a special case of the similarity solutions derived in 
^and corresponds to a fixed point in the differential system for ^5(77), m(ry), /(ry)^ given 

by (16.31) . Thus we find that q — Qq = (qg, fho , /o); is given by 


qo = 


and 


25 - _ 25 

1 ^’ ^ 

In terms of the original variables this corresponds to 


/o = 


25 (25”-1) 

4 ^ S 


Q = 


9t 


M = 


9t2 


and F = 


k^{2S-l)z^ 


(Dl) 


(D2) 


Notably this solution is independent of the source flux of buoyancy. 

We may examine the stability of this fixed point in terms of the similarity variable by 
introducing q = Qo + 9i linearizing about the fixed point q^. This gives 


/ 

3 

4 

0 ^ 


/ -3 
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-1 

4S 

3 

0 

d^i 

3(2S-1) 

2S 

3(2S-1) 

S 

0 

3(2S-1) 

\ 4 

3(2S-1) 

2 

3/ 

V 25 

2(1-45) 

8S 

3 / 


(D3) 


We look for a solution of the form q^ = and deduce that 

(A -h 4S')(6452 a2 - 725'2A - 72SX^ + 126S'A - 1625' - 36A -f 81) = 0. (D 4) 


From this condition, we deduce that there is always a root for which Re (A) > 0 when 
5 > 1, while there are other roots for which Re (A) < 0. Thus the fixed point, q^ is 
linearly unstable and this implies that in terms of the similarity variables, while the 
solution may approach the fixed point, it does not remain close to it asymptotically as 
?7 —^ 00. We illustrate this by plotting q/\fm as a function of the similarity variable q for 
J" = 20, 2 and 0.05 (see figure [161). ^ steady state corresponds to q/y/^ = 1, while for 
the separable solution q/y/^ = 5/9. We note that while the evolution for the strongest 
decrease in buoyancy flux {T = 20) becomes close to qjyf^ = 5/9 at one location during 
its evolution, this separable solution does not in general strongly characterise the entire 
form of the similarity solution. 
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Figure 16. The scaled width of the plume, qjV^, as a function of the similarity variable, rj, 
for T — 0.05 (green line, colour online), J- = 2 (red line, colour online) and = 20 (blue line, 
colour online) and shape factor S — 1.1. 
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